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GREEK  LETTERS 

Yi  =  Activity  coefficient  in  the  adsorbed  phase 
e  =  Adsorption  site  energy  (kJ/mole) 

Up  =  p***  central  moment  of  a  probability  density  function 
71  =  Spreading  pressure 

p  =  Site  matching  parameter  or  covariance  of  the  probability  density  function 

a  =  Square  root  of  the  variance  of  a  probability  density  function 

o*  =  Excess  area  of  mixing  (kg/mole) 

(|)  =  Fugacity  coefficient 

i|f  =  Spreading  pressure  function  7tA,/RT  (moles/kg) 
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=  Area  of  the  Adsorbent 
=  Langmuir  isotherm  fitting  parameter  (1/MPa) 

=  Cumulative  (integral)  energy  distribution 
=  Single  probability  density  function 
=  Joint  probability  density  function 
=  Exponential  curve  fitting  parameters 
=  Freundlich  parameters 
= 

=  Global  amount  adsorbed  (moles/kg) 

=  Amount  adsorbed  in  local  homogeneous  patch  (moles/kg) 
=  Partial  pressure  (MPa) 

=  Saturation  adsorption  value  (moles/1^) 

=  Universal  gas  constant  (kJ/mole»K  or  MPa»L/mole»K) 
=  Constant  separation  factor  isotherm  parameter 
=  Temperature  (K) 

=  Adsorbed  phase  mole  fraction 
=  Vapor  phase  mole  fraction 
=  Total  number  of  components  in  the  mixture 


THEORETICAL  ASPECTS  OF  MULTICOMPONENT 
ADSORPTION  EQUILIBRIA 


INTRODUCTION 

TTie  need  exists  for  a  method  to  accurately  model  mixture  adsorption  equilibria.  This 
need  arises  not  only  from  an  academic  desire  to  understand  the  fundamentals  of  the  adsorption 
process,  but  also  from  the  more  pragmatic  requiremrat  to  provide  the  nucleus  of  mathematical 
models  which  describe  dynamic  adsorption  systems.  In  the  course  of  generating  results  for  a 
dynamic  adsorption  system  using  accepted  mathetical  models,  adsorption  iostherm 
relationships  are  regenerated  hundred  of  thousands  of  times.  This  points  to  the  importance  of 
generating  equilibrium  adsorption  data  with  a  minimum  of  computational  burden  and  a  high 
degree  of  accuracy. 

Several  models  have  been  developed  to  describe  adsorption  systems.  The  simpler 
models  often  are  limited  in  their  applicability  and  accuracy,  while  the  more  complicated 
systems  become  too  cumbersome  or  unsolvable  under  some  conditions.  Several  of  the  models 
are  ultimately  correlative  in  nature.  While  these  models  provide  important  information  in 
understanding  the  mixture  adsorption  process,  their  utility  is  limited  due  to  the  paucity  of,  and 
difficulty  in  obtaining,  mixture  adsorption  data.  As  a  result,  models  that  can  predict  mixture 
adsorption  from  single  component  adsorption  data  are  becoming  more  popular. 

This  study  will  investigate  one  of  the  more  popular  mixture  adsorption  theories,  the 
Adsorbed  Solution  Theory  (AST)  (Myers  and  Prausnitz,  1965),  and  its  ability  to  predict 
mixture  adsorption  from  single  component  data.  In  an  effort  to  isolate  pure  theory  from  the 
effects  and  inherent  errors  associated  with  experimental  data,  this  work  will  focus  mostly  on 
theoretical  manipulations  and  evaluations  of  the  theory.  This  report  details  the  development  of 
a  computationally  fast  adsorption  isotherm  model  that  can  be  applied  to  systems  which 
demonstrate  energetic  heterogeneity.  Since  this  work  focuses  primarily  on  the  AST,  a 
background  section  providing  the  details  of  this  theory  is  included  below. 
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I:  ADSORBED  SOLUTION  THEORY 
BACKGROUND 

The  Adsorbed  Solution  Theory  as  proposed  by  Myers  and  Prausnitz  (1965)  is  based  on 
the  assumption  that  an  adsorbed  mixture  can  be  riq)resCTted  as  a  homogeneous,  well-mixed, 
two-dimensional  phase.  This  assumption  leads  to  a  series  of  thermodynamic  relations  confined 
to  a  two-dimensional  plane.  Standard  three  dimensional  properties  (Pressure  and  Volume)  are 
reduced  to  areal  properties  (Spreading  Pressure  and  Area).  Functions  which  describe 
nonidealities,  such  as  activity  coefficients  and  excess  volume  of  mixing,  also  translate  to  two 
dimensional  analogs.  A  consequence  of  this  approach  is  that  adsorbate/adsorbent  interactions 
in  the  AST  are  completely  separate  from  adsorbate/adsorbate  interactions  which  are  modeled 
as  mixture  nonidealties.  The  theory  uses  the  spreading  pressure  of  the  mixture  as  the  reference 
state.  The  resulting  equations  are  as  follows: 


1  ^  X 

-L  =  j;  i  .  o' 

^TOT  '=1  n° 

(1.1) 

2 

1  =  E  -f, 

(1.2) 

=  pp.p^'iA': 

(1.3) 

P°  '  -^.(O 

(1.4) 

<|r  = 

(1.5) 

n.  -  Xj  ^tot 

(1.6) 

The  function  ^ [(nj®)  is  determined  directly  from  the  single  component  isotherm  for  each 
specie.  The  function  .^(n;®)  is  found  using  the  following  equation: 


2 


d\nn 


(1.7) 


/ 


Or  whra  i|r  is  a  function  of  P  the  9'^  integral  can  take  the  form: 


(1.8) 


The  integrand  of  Equation  1.7  or  1.8  is  obtained  using  the  single  component  isotherm  for  each 
specie  using  the  mixture  spreading  pressure  as  the  standard  state.  A  more  rigorous  derivation 
of  the  thermodynamics  associated  with  the  AST  can  be  found  in  Van  Ness  (1969). 

If  it  is  assumed  that  the  vapor  and  adsorbed  phases  are  ideal  (Ideal  Adsorbed  Solution 
Theory  -  lAST)  o*  becomes  0  and  Yi,  <l>i.  and  <|)®  all  become  1.  The  resulting  5Z+2  unknowns 
in  4Z+2  equations  can  be  determined  if  Z  conditions  are  known  (i.e.  partial  pressures  for  each 
component).  Thus,  the  method  can  predict  multi-component  adsorption  values,  requiring  only 
single  component  isotherm  information  for  each  specie  in  the  mixture. 

In  the  lAST,  the  resultant  equations  are  greatly  simplified.  However,  this 
simplification  restricts  the  application  of  the  theory  to  adsorbates  that  behave  ideally  as  a 
mixture.  In  adsorption  systems  that  deviate  significantly  fi’om  ideality,  a  model  which 
describes  the  adsorbed  phase  activity  coefficients  and  excess  area  of  mixing  must  be  employed. 
Costa  et  al.  (1981)  developed  the  Real  Adsorbed  Solution  Theory  (RAST)  to  model  non-ideal 
adsorption  systems  using  liquid  phase  activity  coefficient  models.  Talu  and  Zwiebel  (1986) 
followed  with  a  proof  that  the  adsorbed  phase  activity  coefficients  must  depend  on  the  system 
spreading  pressure  in  order  to  be  thermodynamically  consistent.  Gamba  et  al.  (1989  and 
1990)  then  modified  the  RAST  to  minimize  the  errors  associated  with  low  pressure  regions  and 
lack  of  thermodynamic  consistency.  However,  use  of  this  modified  RAST  is  limited  to 
applications  requiring  calculations  in  a  small  concentration  range  where  binary  adsorption  data 
already  exist. 

H:  THE  DEVELOPMENT  OF  A  FAST  HETEROGENEOUS  BINARY 
ADSORPTION  ISOTHERM 

BACKGROUND 

An  important  practical  purpose  for  an  adsorption  equilibrium  equation  is  to  aid  in  the 
modeling  of  adsorption-based  separation  systems.  The  models  that  describe  such  systems 
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accurately  are  often  complex  and  typically  involve  several  layers  of  nested  iterations.  At  the 
root  of  all  the  calculations  is  the  equilibrium  relation  which  may  get  called  hundreds  of 
thousands  of  times  in  any  givoi  simulation.  Therefore,  a  simple,  fast  equilibrium  relation  is 
desired  to  reduce  the  computational  burdra  and  create  a  manageable  overall  model. 

In  an  effort  to  address  the  need  for  a  computationally  fast  equilibrium  relation,  O'Brien 
and  Myers  (1985)  developed  a  fast  adsorption  equilibrium  algorithm  (FastlAS)  based  upon  the 
Ideal  Adsorbed  Solution  Theory  (Myers  and  Prausnitz,  1965).  Although  this  model  is  many 
times  ^ter  than  the  analogous  iterative  solution,  it  is  restricted  to  homogeneous  adsorption. 

To  complicate  matters  further,  the  single  component  isotherm  relation  used  in  the  FastlAS,  the 
TALAN  (Taylor  series  expansion  of  the  Langmuir  isotherm)  (O'Brien  and  Myers,  1984),  is  a 
heterogeneous  isotherm  model.  This  application  of  a  heterograeous  single  component 
isotherm  relation  in  a  homograeous  mixture  adsorption  model  seems  somewhat  inappropriate. 
Additionally,  since  the  majority  of  adsorption  systems  of  interest  can  be  categorized  as 
heterogeneous,  a  new,  fast  algorithm  is  needed.  This  section  details  the  development  of  a  fast 
heterogenous  binary  adsorption  equilibrium  relation.  Because  of  the  need  for  speed,  only 
those  solutions  that  produce  a  single  analytical  equation  are  sought. 

Adsorbent  heterogeneity  has  been  historically  modelled  by  describing  the  surface  as 
either  patchwise  or  randomly  heterogeneous  (Jaroniec  and  Madey,  1988,  p.8).  Patchwise 
heterogeneity  assumes  that  like  sites  are  grouped  together  to  form  independent  homogeneous 
patches,  while  random  heterogeneity  assumes  that  the  individual  homogeneous  sites  are 
ungrouped  and  randomly  distributed  across  the  surface.  The  patchwise  model  has  been 
favored  by  researchers  in  the  past  because  the  location  or  relative  placement  of  the  sites  is  not 
needed.  In  the  random  model,  site  layout  is  needed  to  determine  lateral  interactions  between 
nearest  neighbors.  In  the  patchwise  model,  since  the  patches  are  homogenous  and 
independent,  all  nearest  neighbors  are  known  to  be  other  homogeneous  sites.  If,  however, 
lateral  interactions  are  ignored,  the  two  models  become  identical. 

For  a  patchwise  heterogenous  adsorbent,  Jaroniec  and  Rudzinski  (1975)  found  that  the 
general  equation  for  heterogeneous  adsoiption  takes  on  a  multiple  integral  form.  In  a  binary 
system  the  equation  becomes: 


Where  N,  is  the  global  amount  of  i  adsorbed,  e,  is  the  energy  of  adsorption,  n,  is  the  local 
amount  of  i  adsorbed  for  a  patch  characterized  by  the  adsorption  site  energies  e,  and  €3,  and 
g(€,,  €3)  is  the  joint  probability  density  function  for  the  distribution  of  adsorption  energies.  A 
similar  equation  can  be  used  to  describe  the  adsorption  of  component  2.  Since  transposing  the 
component  subscripts  in  the  equations  produces  the  component  2  analog,  the  remainder  of  this 
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work  will  consider  only  the  component  1  equation. 

The  choice  of  the  local  isotherm  function  is  determined  by  the  type  of  adsorption 
(monolayer  or  multilayer)  and  lateral  interactions  that  are  expected  on  the  homogeneous  patch. 
The  majority  of  past  research  has  focused  on  monolayer  adsorption  without  lateral  interactions 
(Jaroniec  and  Madey,  1988,  p.  220).  While  Jaroniec  and  Madey  (1988,  pp.  213-215,  220) 
have  investigated  heterogenous  adsorption  using  local  isotherms  with  lateral  interactions  and 
multi-layer  adsorption,  the  equations  are  too  cumbersome  to  produce  an  analytical  solution.  In 
addition,  as  Jaroniec  and  Madey  (1988,  p.  220)  state,  there  is  no  satisfactory  theory  for 
multilayer  adsorption  that  can  used  as  a  homogeneous  local  isotherm.  For  these  reasons, 
this  work  will  focus  on  monolayer  adsorption  with  no  lateral  interactions. 

The  most  common  adsorption  equation  for  monolayer  adsorption  with  no  lateral 
interactions  is  the  Langmuir  adsorption  isotherm.  If  it  is  assumed  tiiat  the  saturated  molar 
surface  coverage  is  the  same  for  both  components  (this  can  later  be  modified  by  adding  the  Q, 
*  Q2  correction  detailed  below),  the  local  isotherm  then  takes  the  form: 


"1  = 


1  + 


(2.2) 


The  joint  probability  density  function  is  typically  developed  by  combining  the 
individual  probability  density  functions  with  a  site  matching  parameter.  The  common 
statistical  representation  of  tiiis  site  matching  parameter  is  the  covariance  of  the  two 
distributions  p.  p  can  take  on  any  value  from  -1  to  + 1  and  is  a  measure  of  how  the  individual 
distribution  functions  are  related:  p=-l  indicates  perfect  negative  correlation,  p=0  indicates 
random  site  matching,  p =  + 1  indicates  perfect  positive  correlation.  Perfect  positive 
correlation  means  that  the  site  with  the  highest  adsorption  energy  for  component  1  will  have 
the  highest  adsorption  energy  for  component  2  with  each  subsequent  site  behaving  in  the  same 
manner  down  to  the  lowest  energy  site  for  each  component.  Perfect  negative  correlation 
means  that  the  site  with  the  highest  adsorption  energy  for  component  1  will  have  the  lowest 
adsorption  energy  for  component  2  with  each  subsequent  site  behaving  in  the  same  contrary 
manner.  Random  site  matching  means  that  for  each  adsorption  site,  the  adsorption  energy  for 
component  1  has  no  bearing  on,  and  is  completely  uncorrelated  with,  the  adsorption  energy  for 
component  2.  All  other  values  of  p  indicate  varying  degrees  of  positive  or  negative  site 
matching. 

The  majority  of  past  research  into  the  predetermined  selection  of  covariance  has  come 
from  the  two  research  groups  associated  with  Jaroniec  and  Myers.  Valenzuela  et  al.  (1988) 
have  indicated  that,  since  most  adsorption  phenomena  are  controlled  by  dispersion  forces, 
adsorbates  that  are  similar  in  size  and  shape  should  behave  similarly.  Therefore,  most  mixture 
adsorption  systems,  in  the  absence  of  chemical  interactions,  are  perfectly  correlated  (p  =  l). 
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Indeed,  Hoory  and  Prausnitz  (1967)  found  that,  for  mixtures  of  ethane  and  ethylene  on 
charcoal,  a  p=l  value  gave  an  excellent  fit  of  the  data.  Jaroniec  (Jaroniec  and  Madey,  1988, 
pp.  199-202,  203-208)  has  also  investigated  perfectly  correlated  adsorption,  but  limited  their 
efforts  to  adsorbates  which  demonstrate  identical  single  component  energy  distributions,  offset 
only  by  the  differences  in  the  means  of  the  individual  distributions.  The  remainder  of 
Jaroniec's  work  in  this  area  focuses  on  the  p  =0  case. 

What  follows  is  an  evaluation  of  the  possible  solutions  to  Equation  1.  In  an  effort  to  be 
comprehensive,  solutions  are  sought  for  each  type  of  site  matching:  p=0,  p=+l  (p=-l 
producing  mathematically  similar  results  to  p= + 1),  and  also  for  p  as  a  binary  fit  parameter. 

In  addition,  two  methods  of  solution  (direct  integration  and  expansion)  are  investigated.  The 
evaluation  is  separated  by  the  solution  approach  and  is  subdivided  by  the  method  of  site 
matching. 


THEORY 

Direct  Integration  Solution 

A  direct,  analytical  solution  from  integration  is  the  preferred  method  of  solution  for 
Equation  2.1.  Other  methods  such  as  numerical  integration,  mathematical  approximations, 
and  expansions  prior  to  integration  may  introduce  unwanted  errors  into  the  calculation. 
Furthermore,  a  numerical  integration  or  inclusion  of  an  extended  series  defeats  the  purpose  of 
developing  a  fast,  analytical  solution  to  the  binary  heterogeneous  adsorption  equation. 

A  direct  integration  requires  the  selection  of  a  joint  probability  density  function.  The 
simplest  function,  the  bivariate  normal  distribution  used  by  Hoory  and  Prausnitz  (1967),  takes 
the  form: 


with: 


- 1  g  2(1  -  P') 

2710^02^(1  -p2) 


(2.3) 


(2.4) 


p  as  a  Fit  Parameter.  Allowing  p  to  be  a  fit  parameter  in  the  binary  heterogeneous 
adsorption  equation  is  appealing  because  it  removes  any  presuppositions  as  to  the  relative 
interaction  of  each  component  with  the  surface.  It  does,  however,  add  the  complexity  that 
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binary  data  are  required  to  determine  p.  In  the  strictest  srase,  the  model  becomes  a  correlative 
rather  than  a  predictive  solution.  Substituting  Equations  2.2  -  2.4  into  Equation  2.1  produces 
a  double  integral  equation  which  is  horribly  complex.  However,  with  a  little  maiupulation  and 
substitution,  the  timer  (component  1)  int^ml  can  be  reduced  to  the  form: 


/— 

^  F  * 


(At*  *Bx) 


F  +  Ge^ 


-dx 


(2.5) 


with  X  representing  Cj  and  A,  B,  F,  G,  and  D  being  complex  functions  of  Cj,  p,  and  assorted 
constants  from  the  above  equations.  Unfortunately,  no  analytical  solution  to  this  integral  has 
been  found. 


p=0  —  Random  Site  Matching.  If  random  site  matching  is  assumed  (p=0).  Equation  2.3 
takes  the  form: 


2nOjOj 


(2.6) 


Unfortunately,  using  Equations  2.2  and  2.6  to  solve  Equation  2.1  gives  an  equation  of  the 
same  form  as  Equation  2.5.  As  with  the  case  where  p  is  a  variable,  no  solution  has  been 
found. 

p=l  —  Perfect  Positive  Correlation.  If  perfect  positive  correlation  is  assumed  (p= +1),  then 
the  joint  probability  density  function  becomes  the  individual  probability  density  function  of  one 
of  the  components.  Valenzuela  et  al.  (1988)  used  this  information  to  simplify  Equation  1  to  a 
single  integral  form  with  respect  to  the  reference  component  1: 

~  f  (2.7) 


Since  the  energy  distributions  of  the  two  components  are  positively  correlated,  Cj  can  be 
represented  as  a  function  of  the  reference  energy  Cj.  According  to  Valenzuela,er  al.  (1988), 
this  function  (ej*)  takes  the  form: 


^2(^2)  =  ^1(^1)  (2.8) 

The  function  F  represents  the  cumulative  (integral)  energy  distributions  for  each  component. 
Equation  2.8  states  that  the  value  of  can  be  found  by  establishing  the  e,  value  for  a  given 
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site,  finding  the  coiie^nding  area  under  the  component  1  probability  density  Ainction  defined 
by  £„  and  choosing  the  component  2  raeigy  value  which  gives  an  equivalent  area  under  the 
component  2  probability  density  function.  If  the  energy  distributions  are  balanced  around  the 
means  (i.e.  not  skewed),  then  ^uation  2.8  simplifies  to  the  linear  form: 


o,  —  _ 

“(^1  “  €j)  +  €2 


(2.9) 


Jaroniec  and  Madey  (1988,  pp.  201-202,  203-208)  performed  a  similar  simplification,  but 
restricted  the  application  to  identically  shaped  energy  distributions  by  fixing  the  slope  of  the 
linear  relationship  to  one.  Unfortunately,  as  with  the  other  cases  noted  above,  substituting 
Equations  2.2  and  2.9  into  Equation  2.7  produces  the  familiar  integral  form  found  in  Equation 
2.5. 

O'Brien  and  Myers  Expansion  Procedure 

A  second  approach  to  solving  Equation  2.1  follows  the  Taylor  series  expansion 
procedure  O'Brien  and  Myers  (1984)  used  to  develop  the  TAIAN  single  component 
heterogeneous  isotherm.  The  authors  expanded  the  energy  dq>endence  of  the  local  Langmuir 
isotherm  about  the  mean  of  the  energy  distribution.  This  approach  is  appealing  because  it  does 
not  require  an  a  priori  specification  of  a  distribution  function.  In  addition,  even  though  this 
method  uses  an  expansion  to  approximate  a  portion  of  the  equation,  the  penalty  for  truncation 
at  each  level  is  well  defined.  Each  additional  term  represents  an  additional  central  moment  of 
the  probability  density  function.  Therefore,  truncating  the  series  at  any  given  point  results  in 
the  specification  of  the  general  form  of  the  probability  density  function.  For  example, 
truncation  of  the  O'Brien  and  Myers  expansion  after  the  second  term  (as  in  the  TALAN)  fixes 
the  shape  of  the  probability  function  to  the  symmetric,  normal  distribution  form.  Truncation 
after  the  third  term  allows  for  skewness.  Truncation  after  the  fourth  term  allows  for  kurtosis. 

It  should  be  noted,  however,  that  each  additional  term  in  the  series  may  provide  a  superior  fit 
to  the  data,  but  adds  an  additional  fit  parameter.  And,  with  each  additional  fit  parameter,  it 
becomes  more  difficult  to  deterrnine  whether  a  superior  fit  is  provided  because  the  energy 
distribution  is  better  described  by  the  additional  moments,  or  because  more  parameters  can 
make  a  fundamentally  poor  model  seem  better. 

To  apply  the  approach  of  O'Brien  and  Myers  to  Equation  2.1,  their  theory  must  be 
extended  to  account  for  more  than  one  component.  In  the  case  of  a  binary  isotherm,  the  local 
(binary)  isotherm  needs  to  be  expanded  about  the  means  of  the  energy  distributions  of  both 
components  yielding: 
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(2.10) 


••  •  c.  \(c  —  c.  ^ 


^=0  ^*0 


”1  (ei.eaXe,  -  e/(e,  - 
(p+9)! 


where  is  the  multiple  derivative  of  the  local  isotherm  with  respect  to  €,  (p*  derivative) 
and  €2  (q"*  derivative).  If  we  follow  O'Brim  and  Myers,  substitute  Equation  2. 10  into 
Equation  2.1  and  then  interchange  the  order  of  integration: 


.  -  „<PX9). 


Af.  =  EE 

/>=0  ^=0 


n 


1 


'Xe,.€2) 


(p+^')! 


f  f  (€,  -  €i)P(€2  -  €2)«g(€i,C2)rfe,^/ej 


(2.11) 


And,  as  with  direct  integration,  the  possible  solutions  can  be  separated  by  the  method  used  to 
match  the  sites. 

p  as  a  Fit  Parameter.  Using  the  distribution  defined  in  Equation  2.3  and  a  variable  p  value, 
g(€i,  €2)  cannot  be  broten  into  two  indq)endent  fimctions  of  Cj  and  €2.  However,  the 
distribution  can  be  divided  such  that  the  inner  (component  1)  integral  becomes: 

(€2  -  (€,  -  (2.12) 


Since  the  integral  is  with  respect  to  €j,  gX^i,  €2)  becomes  essentially  fi(ej).  Thus,  the  integral 
in  Equation  2. 12  is  the  p*  central  moment  of  fi(€i)  or  Pp.  Equation  2.11  then  becomes: 


Af.  =  EE 

p=0  ^=0 


n 


(PX9)/— — ' 


(^1.€2) 


(p+q)\ 


f  V^p(^2  - 


(2.13) 


Unfortunately,  Pp  is  a  function  of  and  therefore  cannot  be  removed  from  the  integral.  Thus, 
no  solution  to  this  equation  has  been  found. 

p=0  -  Random  Site  Matching.  If  random  site  matching  is  assumed,  the  joint  probability 
density  function  can  be  separated  into  two  independent  functions  as  shown  in  Equation  2.6. 
These  separate  functions  fi(ej)  and  f2(62)  can  be  directly  substituted  for  g(e„e2)  in  Equation 
2.11.  Following  the  work-up  of  Equations  2. 12  and  2. 13,  and  recognizing  the  p*^  and  q*  - 
central  moments  leads  to  the  equation: 
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(2.14) 


p*0q*0  (P^^y 


ptp  and  ft,  arc  respectively  the  p*  central  moment  of  the  component  1  distribution  and  the  q* 
central  moment  of  the  component  2  distribution.  Following  O'Brien  and  Myers  (1984),  = 

1  and  ft]  =  0,  so  that  Equation  2. 14  can  be  written: 


-  -w, 


=  «j(«1.€2)  +  E- 

i»-2 


pi 


-  -(«), 

E- 

«-2 


9! 


EE 

p-2  ,-2 


_(px«)  — .. 

"I 


(p*q)\ 


(2.15) 


The  sums  can  be  extended  knowing  that  ftj/o*  is  the  variance  of  the  distribution,  while  ftj/o’ 
rq)resents  the  skewness  and  fijd*  represents  the  kurtosis. 


Using  Equation  2.2  for  the  local  isotherm  and  following  O'Brien  and  Myers  (1984)  in 
their  development  of  the  TALAN  (truncating  the  sums  after  the  second  terms)  gives: 


(2.16) 


where: 


iL  ii 

$  =  1  +  X  +  ^ 


(2.17) 


It  can  be  seen  that  Equation  2.16  reduces  to  the  TALAN  for  single  component  (P2  =  y  =  0): 


N,  =  0( 


X 

(1-.X) 


2R^T\\+xf 


(2.18) 


Indeed,  fitting  each  component  to  a  TALAN  isotherm  gives  the  parameters  (b,  Q,  e,  and  o) 
for  the  binary  equation.  In  addition,  assuming  a  homogeneous  surface  (Oj  =  Oj  =  0,  Cj  =  Cj, 
and  €2  =  €2)  yields  the  familiar  Langmuir  (Equation  2.2)  form  for  both  binary  and  single 
component  systems. 
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pasi  —  Perfect  Positive  Correlation.  As  shown  above,  in  the  case  of  perfect  positive 
correlation,  Equation  2. 1  amplifies  to  Equation  2.7.  If  Equation  2.9  is  used  to  rq)lace  the 
dependence.  Equation  2.7  can  be  expanded  a  single  time  about  the  mean  of  the  component  1 
raergy  distribution  giving: 


,  £  r  (2.19) 

The  integral  in  Equation  2.19  is  the  p*-central  moment  of  the  component  1  energy  distribution. 
To  evaluate  Equation  2. 19,  Equation  2.2  is  used  as  the  local  isotiiCTm.  Prior  to  differentiation 
in  Equation  2.19,  the  Cj  dependence  is  rq)laced  in  Equation  2.2  using  the  Cj  function  in 
Equation  2.9.  The  sum  is  then  evaluated  at  and  §2-  Truncating  the  sum  after  the  second 
term  gives: 

AT,  .  fi(£  *  ..(1  -  3x  ^  my(m  ,  2(x  *  (2.20) 


where: 

X  =  Pyb^e 


As  with  the  p=0  case,  Equation  2.20  reduces  to  the  TALAN  (Equation  2.18)  for  single 
component  (P2=  y  =  0),  and  to  the  basic  Langmuir  for  a  homogeneous  surface  (Oj  =  o,  =  0, 
e,  =  €i,  and  €2  =  €2). 

Binary  Systems  with  Different  Molar  Monolayer  Coverages  (Q1  *  Q2) 

Equations  2.16  and  2.20  represent  analytical  solutions  for  a  binary  heterogeneous 
system  which  adsorbs  locally  as  a  monolayer  without  lateral  interactions.  These  solutions 
cannot,  however,  be  applied  when  the  molar  monolayer  coverages  for  the  two  components  are 
not  equal.  LeVan  and  Vermeulen  (1981)  addressed  this  problem  with  binary  homogeneous 
adsorption  using  a  Taylor  series  expansion  about  the  mean  of  the  two  different  saturation 
(monolayer)  coverages.  Frey  and  Rodrigues  (1994)  extended  this  solution  to  allow  for  more 
than  two  components.  Since  this  research  is  focusing  on  binary  adsorption,  only  the  LeVan 
and  Vermeulen,  binary  solution  will  be  considered  here.  The  equations  developed  from  the 
LeVan  and  Vermeulen  effort  include  both  a  two-term  and  three-term  Taylor  expansion.  The 


y  = 


RT 


^  =  1  +  X  +  y 


(2.21) 
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primary  equation  for  the  three-term  expansion  takes  the  form: 


n 


* 


(2.22) 


with: 


.  (Cl  -  ft) 

(2,  +e,)x./  x-y  (2.23) 

^  3x^  +  4x  +  xy  -  2y  -  ly^y. 


(x  +  yf 


1  ,y^  +  2v  -  4x  -  x^ 


Q,x  *  1 

(Qx*  02)(x  +  yyx  +  j; 


*  |•)ln(0  -  1) 


(2.24) 


For  the  two  term  expansion,  the  following  simplifications  to  the  above  equations  apply: 


^  ^  Qix  ^ 

X  +  y 


=  0 


(2.25) 


This  analytical  solution  for  different  monolayer  coverages  assumes  ideal  adsorbed  solution 
theory  thermodynamic  interactions.  As  a  result,  this  solution  may  be  used  as  the  local 
isotherm  in  the  proposed  heterogeneous  model.  The  O’Brien  and  Myers  (1984)  expansion, 
detailed  in  equation  2,11,  can  be  applied  using  equation  2.22  as  the  local  isotherm,  and  the 
ultimate  solution  obtained.  Due  to  the  complexity  of  the  resultant  equations,  only  the  two- 
term  Taylor  expansion  form  of  Equation  2.22  will  be  considered. 

p=0  —  Random  Site  Matching.  When  random  site  matching  is  assumed.  Equation  2.22  is 
used  for  the  local  isotherm  in  Equation  2,15.  Truncating  the  expansion  after  the  second  term 
and  using  the  definitions  of  Equations  2.17,  2.23,  2.24,  and  2.25  gives: 
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(2.26) 


+  A 


'U 


2R^T^ 


(o  +  p)  +  Y 


2R^T\x*yf 


where: 


3x^(gi  -0  ^  2x\Q-Q^)  2x\Q-Q,)  3Qx\  2Qx^ 

(*+y)5  (x+yfi  (x+yy^  ^  V  V 


(2.27) 


xyiQ^-Q)  ^  ^\Q-Q2)  ^  ^\Q-Q2)  _Qxy^  2Qjo^^ 
(x+y)^  (x+yf^  (x+y)^  ^  ^ 


(2.28) 


Y 


(o5+o^ln$  +(jco?+yo5)(|--^-)+ 


2„^^^.2„hr 


(x+yf  (x*y)^  ^ 


(2.29) 


p= +1  —  Perfect  Positive  Correlation.  When  perfect  positive  correlation  is  assumed, 
Equation  2.22  is  used  for  the  local  isotherm  in  Equation  2.19.  Truncating  the  expansion  after 
the  second  term  and  using  the  definitions  of  Equations  2.21,  2.23,  2.24,  and  2.25  gives: 


«1  rX 

- [—a  + 

2R^T^  5 


(x+y)^ 


(P-Y)] 


(2.30) 


where: 


5  5  5^ 


(2.31) 


P 


w(/n+l)^ln(5) 


^  (x+otV)+2(x+ot>^)(7w+1)  (x+zny)^ 


(2.32) 


13 


(2.33) 


Y  =  ^  2iii(2x+3iii>>->)ln(g) 

(x+y)5  {x+y) 


'o'  - 

(x^y) 


(2.34) 


W'  =  _  2{x*my)-^, 

{x*y)  (x+y)  ^ 


(2.35) 


It  can  be  seen  that,  with  a  little  mathematical  manipulation,  Equations  2.26  and  2.30  reduce  to 
the  TALAN  (Equation  2.18)  in  the  limit  of  one  component  (P2  =  y  =  0),  to  the  Langmuir 
(Equation  2.2)  when  the  surface  is  homogeneous  (o,  =  Oj  =  0,  Cj  =  e„  and  €2  =  Cj).  and  to 
Equations  2.16  and  2.20  when  the  monolayer  coverages  are  equal  (Qi  =  02  =  Q).  These 
equations  may  be  extended  to  a  higher  number  of  components  or  to  a  greater  level  of 
complexity  in  any  of  the  Taylor  expansions. 

COMPARISONS  AND  DISCUSSION 

Equations  2.16,  2.20,  2.26  and  2.30  rq)resent  a  step  forward  in  the  development  of  a 
computationally  fast  heterogeneous  adsorption  equation.  However,  what  ultimately  determines 
the  utility  of  the  equations  and  underlying  assumptions  is  the  quality  of  the  fit  produced  when 
compared  with  an  exact  solution.  Since  a  universal  exact  solution  to  Equation  1  has  not  been 
found,  the  "exact"  solution  presented  for  each  case  will  be  one  found  when  assuming  specific 
constraints  on  the  system.  What  follows  is  a  discussion  of  how  the  developed  equations 
compare  to  the  "exact"  solutions  calculated  for  each  specific  case. 

Equal  Monolayer  Coverages:  Qj  =  Qj 

p=0  ~  Random  Site  Matching.  Exact  solutions  were  generated  and  compared  to  the  results 
obtained  from  Equation  2. 16  for  a  variety  of  conditions.  Results  are  plotted  in  Figures  2.1- 
2.3.  Results  from  the  Langmuir  isotherm  (Equation  2.2  with  e,  =  e,  and  were  also 

compared  with  the  exact  solutions.  The  exact  solutions  were  generated  by  substituting 
Equations  2.2  and  2.6  into  Equation  2.1  and  integrating  the  results  numerically  using  the  Euler 
method.  It  should  be  noted  that  the  "exact"  solution  is  constrained  by  the  assumptions  inherent 
in  Equations  2.2  and  2.6  (i.e.  local  Langmuir  isotherm  and  Gaussian  shaped  energy 
distribution  curves).  The  Euler  integration  was  centered  on  e  for  each  component  using  200 
integration  steps.  The  selected  model  isotherm  parameters  were  based  on  typical  values  as 
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indicated  by  O’Brien  and  Myers  (1984)  and  are  listed  in  Table  2. 1.  A  copy  of  the  BASIC 
code  writtei  to  perform  the  integration  and  compare  the  results  to  the  calculated  Equation  2.16 
and  Langmuir  results  can  be  found  in  Appendix  A. 

TABLE  2.1:  Variable  values  used  in  numerical  calculations 


Component  1 

Component  2 

T 

298  K 

298  K 

Q 

1.0 

1.0 

€ 

10*RT 

12*RT 

b 

1.0*e-“ 

1.2*e“ 

Figure  2. 1  compares  the  percent  error  for  Equation  2. 16  and  the  Langmuir  isotherm  as  a 
function  of  total  system  pressure  and  breadth  of  the  heterogeneous  energy  distributions  (o/RT). 

In  this  plot,  o/RT  and  the  partial  pressure  for  each  component  are  held  equal  (Yi  =  0.5).  As  the 
plot  shows,  the  deviations  from  the  exact  solution  for  Equation  2.16  are  consistently  smaller  than 
those  of  the  Langmuir  equation.  While  the  Langmuir  errors  get  consistently  worse  with 
increasing  pressure,  the  Equation  2.16  errors  appear  oscillatory.  It  cannot  be  determined  from  the 
plot,  however,  whether  the  Equation  2. 16  series  will  ultimately  converge  at  higher  pressures. 

Figures  2.2  and  2.3  show  the  dramatic  difference  between  the  Equation  2. 16  solution  and 
the  Langmuir  solution.  Figure  2.2  ^ves  the  Equation  2. 16  and  Langmuir  errors  when  the  energy 
distributions  for  both  components  are  held  constant  and  equal  (o/RT  =  0.75)  while  the  relative 
vapor  phase  concentrations  are  varied.  Figure  2.3  compares  Equation  2.16  and  Langmuir  errors 
when  Y,  and  Oj/RT  are  held  constant  at  0.5  and  0.75  respectively,  and  Oj/RT  is  varied.  In  both 
plots,  the  Langmuir  errors  increase  consistently  with  increasing  pressure  while  the  Equation  2. 16 
errors  oscillate  around  0.  As  the  plots  show,  for  any  given  point,  the  errors  obtained  for  the 
Equation  2.16  solution  are  approximately  an  order  of  magnitude  better  than  those  obtained  from 
the  Langmuir  solution. 

p=l  —  Perfect  Positive  Correlation.  Exact  solutions  were  generated  and  compared  to  the 
results  obtdned  from  Equation  2.20  for  a  variety  of  conditions  and  plotted  in  Figures  2.4  -  2.6. 
The  exact  solutions  were  generated  by  substituting  Equations  2.2  and  2.9  into  Equation  2.7  and 
integrating  the  results  numerically  using  the  Euler  method.  As  with  the  p=0  case  above,  the 
"exact"  solution  is  constrained  by  the  assumptions  inherent  to  Equations  2.2  and  2.9  (i.e.  local 
Langmuir  isotherm  and  Gaussian  or  non-skewed  energy  distribution  curves).  The  Euler 
integration  was  centered  on  Cj  using  200  integration  steps.  The  selected  model  isotherm 
parameters  are  the  same  as  those  used  in  the  p=0  calculations  above,  and  can  be  found  in  Table 
2.1.  A  copy  of  the  BASIC  code  written  to  perform  the  integration  and  compare  the  results  to  the 
calculated  Equation  2.20  and  Langmuir  results  can  be  found  in  Appendix  B. 
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FIGURE  2.1:  Percent  Error  for  the  Langmuir  Isotherm  and  Equation 
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FIGURE  2.3:  Percent  Error  for  the  Langmuir  Isotherm  and  Equation 
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FIGURE  2.4:  Percent  Error  for  the  Langmuir  Isotherm  and  Equation  2.20 
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FIGURE  2.5:  Percent  Error  for  the  Langmuir  Isotherm  and  Equation  2.20 
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FIGURE  2.6:  Percent  Error  for  the  Langmuir  Isotherm  and  Equation  2.20 


Figure  2.4  compares  the  percent  error  for  Equation  2.20  and  the  Langmuir  isotherm  as  a 
function  of  total  system  pressure  and  breadth  of  the  heterogeneous  energy  distributions  (o/RT). 

In  this  plot,  o/RT  and  the  partial  pressure  for  each  component  are  held  equal  (Y,  =  0.5).  As  the 
plot  shows,  the  deviations  from  the  exact  solution  for  ^nation  2.20  are  consistently  smaller  than 
those  of  the  Langmuir  equation.  As  would  be  e)q)ected  for  an  expansion  solution,  the  deviations 
from  the  exact  values  are  osdllatory  with  decreasing  absolute  magnitudes  as  the  breadth  of  the 
distribution  gets  narrower.  Although  the  error  gets  as  high  as  8%  for  the  o/RT  =  1.25  solution, 
the  results  still  converge  with  increasing  pressure  values. 

Figure  2.5  shows  how  Equation  2.20  and  the  Langmuir  isotherm  perform  when  the  energy 
distributions  for  both  components  are  held  constant  and  equal  (o/RT  =  0.75)  while  the  relative 
vapor  phase  concentrations  are  varied.  Once  again.  Equation  2.20  ^ves  consistently  better  results 
than  the  Langmuir  isotherm. 

The  final  plot  in  this  group.  Figure  2.6,  provides  the  most  insight  into  the  quality  of  the 
Equation  2.20  results  as  compared  to  those  of  Ae  Langmuir  equation.  In  this  plot,  Yi  and  Oi/RT 
are  held  constant  at  0.5  and  0.75  respectively,  while  Oj/RT  is  varied.  As  Oj/RT  drops  from  0.75 
to  0.25,  the  Equation  2.20  error  shifts  upward  in  the  lower  pressure  region,  but  converges  more 
rapidly  as  the  pressure  is  inCTeased.  The  Langmuir  solution,  however,  is  so  poor  that  the  Oj/RT  = 
0.25  curve  does  not  fit  on  the  plot. 

Unequal  Monolayer  Coverages:  Qi  ^  Qj 

p=0  —  Random  Site  Matching.  Figures  2.7-2.12  are  plots  comparing  Equation  2.26  results  to 
those  obtained  from  the  homogeneous  Langmuir  solutions  of  LeVan  and  Vermeulen  (Equation 
2.22).  The  results  from  each  of  these  equations  are  plotted  as  a  percentage  of  error  calculated 
from  an  "exact"  solution.  The  exact  solutions  were  generated  by  substituting  Equation  2.6  and  an 
ideal  adsorbed  solution  theory,  iterative  mixture  solution  based  on  Equation  2.2,  into  Equation 
2. 1  and  integrating  the  results  numerically  using  the  Euler  method.  As  with  the  Qi  =  Qj  case 
detailed  previously,  the  "exact"  solutions  are  subject  to  the  constraints  inherent  in  Equations  2.2 
and  2.6  (i.e.  local  Langmuir  isotherm  and  Gaussian  shaped  energy  distribution  curves).  The  Euler 
integration  was  centered  on  e  for  each  component  using  200  integration  steps.  The  selected 
model  isotherm  parameters  are  listed  in  Table  2. 1 .  A  copy  of  the  BASIC  code  written  to  perform 
the  integration,  and  compare  the  results  to  the  calculated  Equation  2.22  and  Equation  2.26  results 
can  be  found  in  Appendix  C. 

The  plots  can  be  divided  into  two  series;  Figures  2.7  -  2.9  using  Qj  =  0.9  and  Qj  =  1.1, 
and  Figures  2. 10  -  2. 12  using  Qj  =  0.8  and  Q2  =  1.2.  The  latter  of  the  two  cases  should  represent 
a  fairly  extreme  spread  in  saturation  values,  with  Qj  being  50%  larger  than  Qj.  As  the  plots 
indicate,  as  values  for  site  energy  distribution  breadth,  vapor  phase  mole  fraction,  and  site  energy 
distribution  ratios  are  changed. 
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FIGURE  2.7:  Percent  Error  Equations  2.22  and  2.26 
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FIGURE  2.8:  Percent  Error  Equations  2.22  and  2.26 
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FIGURE  2.9:  Percent  Error  Equations  2.22  and  2.26 


FIGURE  2.10:  Percent  Error  Equations 


(:^U0OJ0ci)  jojjg 


27 


FIGURE  2.11:  Percent  Error  Equations  2.22  and  2.26 
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FIGURE  2.12:  Percent  Error  Equations 


definite  error  patterns  develop  ^ch  are  consistent  for  both  the  Equation  2.22  and  Equation  2.26 
solutions.  The  major  difference  between  the  two  solutions  is  that  the  Equation  2.26  solution 
consistently  provides  a  corrective  offset.  This  offset  makes  the  Equation  2.26  solution 
significantly  superior  (as  much  as  an  order  of  magnitude)  to  the  Equation  2.22  solution.  In 
addition,  as  would  be  expected  for  a  Taylor  series  expansion  solution,  as  the  spread  in  Q  values 
\wdens,  the  solutions  for  both  equations  become  less  accurate  and  more  divergent.  It  is  expected 
that  this  trend  would  continue  if  mixtures  demonstrating  larger  differences  in  Q  values  were 
investigated. 

p=l  -  Perfect  Positive  Correlation.  Figures  2.13-2.18  are  plots  comparing  Equation  2.30 
results  to  those  obtained  fi'om  the  homogeneous  Langmuir  solutions  of  LeVan  and  Vermeulen 
(Equation  2.22).  The  results  from  each  of  these  equations  are  plotted  as  a  percentage  of  error 
calculated  fi'om  an  "exact"  solution.  The  exact  solutions  were  generated  by  substituting 
Equations  2.9  and  2.22  into  Equation  2.7  and  integrating  the  results  numerically  using  the  Euler 
method.  As  det^led  previously,  the  "exact"  solutions  are  subject  to  the  constraints  inherent  in 
Equations  2.7  and  2.22  (i.e.  local  Langmuir  isotherm  and  Gaussian  shaped  energy  distribution 
curves).  The  Euler  integration  was  centered  on  i  for  each  component  using  200  integration  steps. 
The  selected  model  isotherm  parameters  are  listed  in  Table  2.1.  A  copy  of  the  BASIC  code 
written  to  perform  the  integration,  and  compare  the  results  to  the  calculated  Equation  2.22  and 
Equation  2.30  results  can  be  found  in  Appendix  D. 

The  plots  can  be  divided  into  two  series:  Figures  2. 13  -  2. 15  using  =  0.9  and  Qj  =  1. 1, 
and  Figures  2. 16  -  2.18  using  Q,  =  0.8  and  Qj  =1.2.  As  with  the  p  =  0  case  outlined  above,  as 
values  for  site  energy  distribution  breadth,  vapor  phase  mole  fi'action,  and  site  energy  distribution 
ratios  are  changed,  definite  error  patterns  develop  which  are  consistent  for  both  the  Equation  2.22 
and  Equation  2.30  solutions.  However,  unlike  the  p  =  0  case,  the  corrective  offset  provided  by 
the  heterogeneity  terms  inherent  in  the  Equation  2.30  solution  does  not  always  provide  a  superior 
solution.  This  offset  consistently  shifts  the  homogeneous  (or  Equation  2.22)  solutions  to  the 
more  positive  errors.  Since  the  error  function  was  defined  as:  moles  adsorbed  actual  -  moles 
adsorbed  theory,  positive  error  readings  indicate  that  the  theory  is  under-predicting  the  amount 
adsorbed.  As  the  plots  show.  Equation  2.30  gives  superior  results  at  the  lower  total  pressure 
(Ptot)  values.  As  Pjot  increases,  both  solutions  drift  toward  more  positive  error  readings.  Thus, 
as  PjoT  is  increased  past  a  threshold  value,  the  Equation  2.30  solution  becomes  inferior  to  that  of 
Equation  2.22. 


CONCLUSIONS 

As  the  various  plots  indicate,  a  non-iterative,  heterogeneous  solution  has  been  developed 
that  provides  predictive  results  superior  to  those  of  analogous  homogeneous  solutions. 
Variations  in  the  model's  constituent  parameters  indicate  that  the  results  hold  for  a  substantial 
range  of  conditions.  The  one  (general)  case  where  the  new  heterogeneous  solution  did  not  out¬ 
perform  the  homogeneous  solution  is  rendered  insignificant  when  the  following  is 
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FIGURE  2.13:  Percent  Error  Equations  2.22  and  2.30 


FIGURE  2.14:  Percent  Error  Equations 
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FIGURE  2.15:  Percent  Error  Equations  2.22  and  2.30 
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FIGURE  2.16:  Percent  Error  Equations  2.22  and  2.30 


FIGURE  2.17:  Percent  Error 
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FIGURE  2.18:  Percent  Error  Equations 


considered:  1)  the  error  for  the  fast  heterogeneous  model,  within  the  range  of  parameters 
investigated,  was  not  appreciably  higher  than  that  obtained  from  the  homogeneous  equation;  2) 
the  parametric  conditions  where  this  occurred  w^  limited  to  higher  total  pressures  and 
divergent  Q  values;  and  3)  the  errors  for  these  specific  cases  were  very  low  relative  to  the 
homogeneous  solution  errors  at  lower  total  pressures.  As  a  result,  when  all  conditions  are 
considered  (i.e.  the  fiill  range  of  pressure  values),  the  new  heterogeneous  model  provides  the 
most  accurate  solution. 

An  important  caveat  to  this  analysis  is  that  the  "exact"  solutions  which  are  used  as  error 
standards  have  inherit  assumptions  that  may  have  had  impact  on  the  results.  In  addition,  the 
"exact"  solution  is  a  theoretical  one  that  is  based  fundamentally  on  a  heterogeneous  adsorption 
system.  Ultimately,  the  utility  of  a  developed  solution  is  dependent  on  its  ability  to  fit  and 
predict  actual  mixture  data. 

An  indication  of  the  importance  of  developing  fast  solutions  can  be  found  in  the  time  it 
took  to  generate  the  "exact"  solutions  in  the  analysis.  On  a  relatively  fast  PC  (80486  running 
at  50  MHz),  the  compiled  code  to  generate  100  data  points  for  each  plot  required  between  10 
and  200  minutes  to  run.  This  computational  burden  can  be  compared  to  the  immeasurably 
short  time  (milliseconds  or  less)  it  takes  to  execute  the  single  calculation  associated  with  the 
fast  heterogeneous  solution.  Considering  that,  as  stated  above,  any  isotherm  subroutine  in  a 
dynamic  model  will  be  called  hundreds  of  thousands  of  times,  it  b^mes  apparent  that  a  fast 
model  of  the  type  developed  here  is  desirable. 
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APPENDIX  A 

BASIC  code  for  Qi  =  Q2  and  p  =  0 
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10  DEFDBL  A-H,  L-Z 
20  DEFINT  I-K 

30  PI  *  3.1415926#:  R  ■  8.312440000000001D-03:  '  KJ/MOLE  K 
40  RT  =  R  *  298 
50  DIMFX(502),  GX(502) 

70  INC  «  200 
80  GTOTl  =  0 

90  INPUT  "ENTER  OUTPUT  FILE  NAME - >  ",  FIL$ 

95  OPEN  FIL$  FOR  OUTPUT  AS  #1 

100  EBl  *  10  *  RT:  EB2  «  12  *  RT 

110  B1  =  EXP(-IO):  B2  *  1.2  *  EXP(-IO) 

115  INPUT  "ENTER  SIOCVl/RT  — >  ",  SRTl 

117  INPUT  "ENTER  SIGMA2/RT  — >  ",  SRT2 

118  SI  =  SRTl  *  RT:  S2  =  SRT2  *  RT 
120  M  *  {S2  /  SI)  2 

130  INPUT  "ENTER  GAS  PHASE  MOLE  FRACTION  Yl - >  ",  Y1 

200  FOR  PTOT  =  0  TO  1!  STEP  .01 

220  Q  =  1 

230  PI  *  Yl  *  PTOT 

240  P2  =  PTOT  -  PI 

260  FX{0)  =  0:  GX(0)  =  0 

290  GCN  =  0 

295  HI  =  2  *  EBl  /  INC:  H2  =  2  *  EB2  /  INC 
300  FOR  I  =  1  TO  INC 
310  El  =  HI  *  I 

320  XS  =  PI  *  B1  *  EXP (El  /  RT) 

325  PSll  =  (El  -  EBl)  /  SI 

340  G1  =  EXP(-(PSI1  ^  2)  /  2)  /  (SI  *  (2  *  PI)  .5) 

345  FCN  =  0 

350  FOR  J  =  1  TO  INC 

360  E2  =  H2  *  J 

370  YS  =  P2  *  B2  *  EXP(E2  /  RT) 

390  PSI2  =  (E2  -  EB2)  /  S2 

500  N1  =  Q  *  XS  /  (1  +  XS  +  YS) 

520  G2  =  EXP(-(PSI2  ''  2)  /  2)  /  (S2  *  (2  *  PI)  .5) 

530  IF  (PTOT=0)  AND  (1=1)  THEN  GTOT2  =  GTOT2  +  G2  *  H2:  GOTO  600 

550  FX(J)  =  N1  *  G2  /  GTOT2 

560  IF  (J  =  INC)  THEN  FCN  =  FCN  +  H2  *  FX(J)  /  3:  GOTO  600 

580  FCN  =  FCN  +  H2  *  ( (2  +  2  *  (J  MOD  2))  *  FX(J))  /  3 

600  NEXT  J 

630  IF  (PTOT  =  0)  THEN  GTOTl  =  GTOTl  +  G1  *  HI:  GOTO  790 

650  GX(I)  =  FCN  *  G1  /  GTOTl 

660  IF  (I  =  INC)  THEN  GCN  =  GCN  +  HI  *  GX(I)  /  3:  GOTO  790 

680  GCN  =  GCN  +  HI  *  ( (2  +  2  *  (I  MOD  2))  *  GX(I))  /  3 

790  NEXT  I 

800  X  =  PI  *  B1  *  EXP (EBl  /  RT) 

810  Y  =  P2  *  B2  *  EXP(EB2  /  RT) 

820  XI  =  1  +  X  +  Y 

850  HH  =  X  *  (1+Y)  *  (1-X+Y)  *S1^2+X*Y*  (Y-X-1)  *  S2  ^,2 

880  HXl  =  Q  *  (X  /  XI  +  HH  /  (2  *  XI  ^  3  *  RT  ^  2) ) 

890  LANG  =  Q  *  X  /  XI 
895  IF  (FCN  =  0)  THEN  GOTO  920 
900  DEL  =  (GCN  -  HXl)  *  100  /  FCN 

910  DEL2  =  (GCN  -  LANG)  *  100  /  FCN 

920  PRINT  USING  "###.####  ";PTOT;GCN;HXl;DEL;LANG;DEL2;GTOTl;GTOT2 

930  PRINT  #1,  USING  "###.#######  PTOT;  DEL;  DEL2 

940  NEXT  PTOT 

945  CLOSE  #1 

950  END 
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APPENDIX  B 

BASIC  code  for  Qi  =  Qj  and  p  =  1 
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10  DEFDBL  A-H,  L-Z 
20  DEFINT  I-K 

30  PI  «  3.1415926#:  R  «  8.312440000000001D-03:  '  KJ/MOLE  K 
40  RT  -  R  *  298 
50  DIM  FX(5002) 

70  INC  »  1000 
80  GTOT  =  0 

90  INPUT  "ENTER  OUTPUT  FILE  NAME  - >  ",  FIL$ 

95  OPEN  FIL$  FOR  OUTPUT  AS  #1 

100  EBl  =  10  *  RT:  EB2  «  12  *  RT 

110  B1  =  EXP(-IO):  B2  =  1.2  *  EXP(-IO) 

115  INPUT  "ENTER  SIOIAI/RT  — >  ",  SRTl 

117  INPUT  "ENTER  SIGMA2/RT  — >  ",  SRT2 

118  SI  =  SRTl  *  RT:  S2  =  SRT2  *  RT 
120  M  =  (S2  /  SI)  2 

130  INPUT  "ENTER  GAS  PHASE  MOLE  FRACTION  Y1 - >  ",  Y1 

200  FOR  PTOT  *  0  TO  1!  STEP  .01 

220  Q  =  1 

230  PI  =  Y1  *  PTOT 

240  P2  =  PTOT  -  PI 

260  FX(0)  =  0 

290  FCN  =  0 

295  H  =  2  *  EBl  /  INC 

300  FOR  I  =  1  TO  INC 

320  El  =  2  *  EBl  *  I  /  INC 

340  E2  =  M  *  (El  -  EBl)  +  EB2 

345  IF  (E2  <  0)  THEN  PRINT  "E2  <  0  ERROR":  END 
350  XS  =  PI  *  B1  *  EXP (El  /  RT) 

360  YS  =  P2  *  B2  "  EXP(E2  /  RT) 

380  PSIl  =  (El  -  EBl)  /  SI 

390  PSI2  =  (E2  -  EB2)  /  S2 

500  N1  =  Q  ♦  XS  /  (1  +  XS  +  YS) 

510  G1  =  EXP(-(PSI1  ^2)  /  2)  /  (SI  *  (2  *  PI)  ^  .5) 

520  G2  =  EXP(-(PSI2  ^  2)  /  2)  /  (S2  *  (2  *  PI)  ''  .5) 

530  IF  (PTOT  =  0)  THEN  GTOT  =  GTOT  +  G1  *  H:  GOTO  700 

550  FX(I)  =  N1  *  G1  /  GTOT 

560  IF  (I  =  INC)  THEN  FCN  =  FCN  +  H  *  FX(I)  /  3:  GOTO  700 

580  FCN  =  FCN  +  H  *  (  (2  +  2  *  (I  MOD  2))  *  FX(I))  /  3 

700  NEXT  I 

800  X  =  PI  *  B1  *  EXP (EBl  /  RT) 

810  Y  =  P2  *  B2  *  EXP(EB2  /  RT) 

820  XI  =  1  +  X  +  Y 

850  HH  =  1  -  (3*X  +  M*Y*(M+2))  /  XI  +  (2*(X+M*Y)  ^2)  /  XI  ^  2 
880  HXl  =  Q  *  (X  /  XI  +  (X  *  HH  *  SI  ^  2)  /  (2  *  XI  *  RT  ''  2)  ) 

890  LANG  =  Q  *  X  /  XI 
895  IF  (FCN  =  0)  THEN  GOTO  920 
900  DEL  =  (FCN  -  HXl)  *  100  /  FCN 

910  DEL2  =  (FCN  -  LANG)  *  lOO  /  FCN 

920  PRINT  USING  "###.#####  PTOT;  FCN;  HXl;  DEL;  LANG;  DEL2;.  GTOT 

930  PRINT  #1,  USING  "###.#######  ";  PTOT;  DEL;  DEL2 

940  NEXT  PTOT 

945  CLOSE  #1 

950  END 
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BASIC  code  for  Qi  Q2  and  p 
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10  DEFDBL  A-H,  K-Z 
20  DEFINT  I-J 
25  DX  «=  .0000000001# 

30  PI  •=  3.1415926#:  R  =  8.312440000000001D-03:  '  KJ/MOLE  K 
40  RT  ■  R  *  298 
50  DIM  FX(202),  GX(202) 

70  INC  «  200 

80  GTOTl  =  0:  GTOT2  =  0 

90  INPUT  "ENTER  OUTPUT  FILE  NAME  - >  ",  FIL$ 

100  EBl  =  10  *  RT:  EB2  *=  12  *  RT 
110  B1  «»  EXP(-IO):  B2  -  1.2  *  EXP(-IO) 

112  INPUT  "ENTER  SIGMAl/RT  - >  ",  SRTl 

113  INPUT  "ENTER  SICT1A2/RT  - >  ",  SRT2 

115  INPUT  "ENTER  Q1  — >  ",  Q1 

117  INPUT  "ENTER  Q2  — >  ",  Q2 

118  SI  »  SRTl  *  RT:  S2  =  SRT2  *  RT 
120  M  =  (S2  /  SI)  2 

130  INPUT  "ENTER  GAS  PHASE  MOLE  FRACTION  Yl  - >  ",  Y1 

200  FOR  PTOT  =  0  TO  1!  STEP  .01 
210  OPEN  FIL$  FOR  APPEND  AS  #1 
230  PI  =  Yl  *  PTOT 
240  P2  =  PTOT  -  PI 
260  FX(0)  =  0:  GX(0)  =  0 

290  GCN  =  0 

295  HI  =  2  *  EBl  /  INC:  H2  =  2  *  EB2  /  INC 
300  FOR  I  =  1  TO  INC 
310  El  =  HI  *  I 

320  XS  =  PI  *  B1  *  EXP (El  /  RT) 

325  PSIl  =  (El  -  EBl)  /  SI 

340  G1  =  EXP  (-(PSIl  ^  2)  /  2)  /  (SI  *  (2  *  PI)  .5) 

345  FCN  =  0 

350  FOR  J  =  1  TO  INC 

360  E2  =  H2  *  J 

370  YS  =  P2  *  B2  *  EXP(E2  /  RT) 

390  PSI2  =  {E2  -  EB2)  /  S2 

500  GOSUB  2000 

520  G2  =  EXP(-(PSI2  ^  2)  /  2)  /  (S2  *  (2  *  PI)  .5) 

530  IF  (PTOT  =  0)  AND  (I  =  1)  THEN  GTOT2  =  GTOT2  +  G2  *  H2 :  GOTO  600 

550  FX(J)  =  N1  *  G2  /  GTOT2 

560  IF  (J  =  INC)  THEN  FCN  =  FCN  +  H2  *  FX(J)  /  3:  GOTO  600 

580  FCN  =  FCN  +  H2  *  ( (2  +  2  *  (J  MOD  2))  *  FX(J))  /  3 

600  NEXT  J 

630  IF  (PTOT  =  0)  THEN  GTOTl  =  GTOTl  +  G1  *  HI:  GOTO  790 

650  GX(I)  =  FCN  *  G1  /  GTOTl 

660  IF  (I  =  INC)  THEN  GCN  =  GCN  +  HI  *  GX(I)  /  3:  GOTO  790 

680  GCN  =  GCN  +  HI  *  ( (2  +  2  *  (I  MOD  2))  *  GX(I))  /  3 

790  NEXT  I 

795  IF  (PTOT  =  0)  THEN  HXl  =  0:  LANG  =  0:  PRINT  GTOTl,  GTOT2:  GOTO  940 
800  X  =  PI  *  B1  *  EXP (EBl  /  RT) 

810  Y  =  P2  *  B2  *  EXP(EB2  /  RT) 

820  XI  =  1  +  X  +  Y 

830  QB  =  (Q1  *  X  +  Q2  *  Y)  /  (X  +  Y) 

835  DL2  =  (Q1  -  Q2)  *  X  *  Y  *  LOG(XI)  /  (X  +  Y)  ''  2 

855  A  =  3  *  X  ^  2  *  (Q1  -  QB)  /  (  (X  +  Y)  *  XI)  +  2  *  X  3  *  (QB  -  Ql) 

/  ((X  +  Y)  2  *  XI)  +  2  *  X  ^  3  *  (QB  -  Ql)  /  (  (X  +  Y)  *  XI  2) 

+  QB*X/XI-3*QB*X"'2/XI^2  +  2*QB*X^3/XI^3 
860  B  =  X  *  Y  *  (Q2  -  QB)  /  (  (X  +  Y)  *  XI)  +  2  *  X  *  Y  ^  2  *  (QB  -  Q2) 

/  (  (X  +  Y)  ^  2  *  XI)  +  2  *  X  *  Y  ^  2  *  (QB  -  Q2)  /  (  (X  +  Y)  *  XI 

^2)  -QB*X*Y/XI''2  +  2*QB*X*Y''2/XI^3 
865  G  =  (SI  ^  2  +  S2  2)  *  LOG(XI)  +  (X  *  SI  ^  2  +  Y  *  S2  ^  2 )  *  (3  / 
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XI  -  6  *  LOG(XI)  /  (X  +  Y))  +  (X  2  *  SI  2  +  Y  2  *  S2  2) 
(6  *  LOG(XI)  /  (X  +  Y)  ^  2  -  4  /  ( (X  +  Y)  *  XI)  -  1  /  XI  ^  2) 


* 


870  • 

875  HH  «  (Q1  -  Q2)  *  X  *  Y  /  (2  *  RT  ^  2  *  (X  +  Y)  '‘2) 

880  HXl  =  QB  *  X  /  XI  +  DL2  +  SI  2  *  (A  +  B)  /  (2  *  RT  ^  2)  +  G  *  HH 

890  LANGL  =  QB  *  X  /  XI  +  DL2 

895  IF  (GCH  «  0)  THEN  GOTO  920 

900  DEL  =  (GCN  -  HXl)  *  100  /  GCN 

910  DEL2  «  (GCN  -  LANGL)  *  100  /  GCN 

915  DEL3  =  (GCN  -  {  (Q1  +  (J2)  /  2)  *  X  /  XI)  *  100  /  GCN 

920  PRINT  USING  "###.#####  PTOT;  GCN;  HXl;  DEL;  LANGL;  DEL2;  DEL3 

930  PRINT  #1,  USING  "###.#######  ";  PTOT;  DEL;  DEL2;  DEL3 

940  CLOSE  #1 

945  NEXT  PTOT 

950  END 

960  ' 

970  ' 

2000  '  LANraiUIR  AST  ALGORITHM 
2100  • 

2110  • 

2115  IF  (PTOT  =  0)  THEN  N1  =  0:  RETURN 
2120  K1  =  B1  *  EXP (El  /  RT) 

2130  K2  *  B2  *  EXP(E2  /  RT) 

2200  IC  =  0 

2210  PSI  =  .0000001 

2220  FOR  IJ  =  0  TO  1  STEP  1 

2230  PSIG  =  PSI  *  (1  +  DX  *  IJ) 

2240  PIO  =  (EXP(PSIG  /  Ql)  -  1)  /  K1 

2250  P20  =  (EXP(PSIG  /  Q2)  -  1)  /  K2 

2260  XI  =  PI  /  PIO 

2270  X2  =  P2  /  P20 

2280  IF  (IJ  =  0)  THEN  ERl  =  XI  +  X2  -  1  ELSE  DERI  =  XI  +  X2  -  1 
2300  NEXT  IJ 

2310  IF  (ABS(ERl)  <  DX  *  1000)  THEN  GOTO  2500 
2320  IC  =  IC  +  1 

2330  IF  (IC  >  50)  THEN  PRINT  "NON  CONVERGE":  END 
2340  DDERl  =  (DERI  -  ERl)  /  (PSI  *  DX) 

2350  PSI  =  PSI  -  ERl  /  DDERl 
2400  GOTO  2220 
2420  ' 

2450  • 

2470  • 

2500  NIO  =  Ql  *  K1  *  PIO  /  (1  +  K1  *  PIO) 

2510  N20  =  Q2  *  K2  *  P20  /  (1  +  K2  *  P20) 

2520  NTOT  =  1  /  (XI  /  NIO  +  X2  /  N20) 

2530  N1  =  NTOT  *  XI 
2550  • 

2600  RETURN 
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APPENDKD 

BASIC  code  for  Qi  *  Qj  and  p 
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10  DEFDBL  A-H,  K-Z 
20  DEFINT  I-J 
25  DX  =  .0000000001# 

30  PI  «  3.1415926#:  R  «  8.312440000000001D-03:  •  KJ/MOLE  K 
40  RT  =  R  *  298 
50  DIM  FX(500) 

70  INC  «=  200 
80  GTOT  »  0 

90  INPUT  "ENTER  OUTPUT  FILE  NAME - >  ",  FIL$ 

95  OPEN  FIL$  FOR  OUTPUT  AS  #1 

100  EBl  «  10  *  RT:  EB2  *  12  *  RT 

110  B1  =  EXP(-IO):  B2  =  1.2  *  EXP(-IO) 

112  INPUT  "ENTER  SIQIAI/RT  - >  ",  SRTl 

113  INPUT  "ENTER  SIOIM/RT - >  ",  SRT2 

115  INPUT  "ENTER  Q1  — >  ",  Q1 

117  INPUT  "ENTER  Q2  — >  ",  Q2 

118  SI  =  SRTl  *  RT:  S2  =  SRT2  *  RT 
120  M  =  (S2  /  SI)  2 

130  INPUT  "ENTER  GM  PHASE  MOLE  FRACTION  Yl - >  ",  Y1 

200  FOR  PTOT  =  0  TO  1!  STEP  .01 

230  PI  =  Yl  *  PTOT 

240  P2  =  PTOT  -  PI 

260  FX{0)  =  0 

290  FCN  =  0 

295  H  =  2  *  EBl  /INC 

300  FOR  I  =  1  TO  INC 

320  El  =  H  *  I 

340  E2  =  M  *  (El  -  EBl)  +  EB2 

345  IF  (E2  <  0)  THEN  PRINT  "E2  <  0  ERROR":  END 

350  XS  =  PI  *  B1  *  EXP (El  /  RT) 

360  YS  =  P2  *  B2  *  EXP(E2  /  RT) 

380  PSll  =  (El  -  EBl)  /  SI 

390  PS12  =  (E2  -  EB2)  /  S2 

500  GOSUB  2000 

510  G1  =  EXP(-(PSI1  ''  2)  /  2)  /  (SI  *  (2  *  PI)  .5) 

520  G2  =  EXP(-(PSI2  ''  2)  (  2)  /  (S2  *  (2  *  PI)  ''  .5) 

530  IF  (PTOT  =  0)  THEN  GTOT  =  GTOT  +  G1  *  H:  GOTO  700 

550  FX(I)  =  N1  ♦  G1  /  GTOT 

560  IF  (I  =  INC)  THEN  FCN  =  FCN  +  H  *  FX(I)  /  3;  GOTO  700 

580  FCN  =  FCN  +  H  *  ( (2  +  2  *  (I  MOD  2))  *  FX(I))  /  3 

700  NEXT  I 

750  IF  (PTOT  =  0)  THEN  HXl  =  0:  LANG  =  0:  GOTO  920 
800  X  =  PI  *  B1  *  EXP (EBl  /  RT) 

810  Y  =  P2  *  B2  *  EXP(EB2  /  RT) 

820  XI  =  1  +  X  +  Y 

830  QB  =  (Q1  *  X  +  Q2  *  Y)  /  (X  +  Y) 

835  DL2  =  (Q1  -  Q2)  *  X  *  Y  *  LOG(XI)  /  (X  +  Y)  ‘^2 

840  QBP  =  ( (Q1  *  X  +  Q2  *  M  *  Y)  -  QB  *  (X  +  M  *  Y) )  /  (X  +  Y) 

850  QBPP  =  (  (Q1*X+Q2*Y*M^2)  -  QB*  (X+Y*M'"2)  -  2*  (X+M*Y)  *QBP)  /  (X  +  Y) 

855  A  =  (QB  +  QBP)  *  (1  -  2  *  {X  +  M  *  Y)  /  XI)  +  (QBP  +  QBPP)  -  QB  * 

(X  +  Y  *  M  2)  /  XI  +  (2  *  QB  *  (X  +  M  *  Y)  ^2)  /  (XI  ^  2) 

860  B  =  M  *  ((M  +  1)  ^2)  *  LOG{XI)  +  ((X  +  Y*M''2)  +  2*  (X+M* 
Y)  *  (M  +  1)  )  /  XI 

865  G  =  (  (X  +  M  *  Y)  ^2)  /  (XI  2)  +  (4  *  (X  +  M  *  Y)  2)  /  (  (X  + 

Y)  *  XI)  +  2*M*  (2*X  +  3*M*Y-Y)  *  LOG(XI)  /  (X  +  Y) 

870 

875  HH  =  X  *  A  /  XI  +  (Q1  -  Q2)  *  X  *  Y  *  (B  -  G)  /  (  (X  +  Y)  ^2) 

880  HXl  =  QB  *  X  /  XI  +  DL2  +  SI  ^  2  *  HH  /  (2  *  RT  ^  2 ) 

890  LANGL  =  QB  *  X  /  XI  +  DL2 

895  IF  (FCN  =  0)  THEN  GOTO  920 
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900  DEL  =  (FCN  -  HXl)  *  100  /  FCN 
910  DEL2  *  (FCN  -  LANGL)  *  100  /  FCN 

915  DEL3  -  (FCN  -  ( (Q1  +  Q2)  /  2)  *  X  /  XI)  *  100  /  FCN 

920  PRINT  USING  "###.#####  ";PTOT; FCN; HXl; DEL; LANGL; DEL2;GTOT;DEL3 

930  PRINT  #1,  USING  "###.#######  ";  PTOT;  DEL;  DEL2;  DEL3 

940  NEXT  PTOT 

945  CLOSE  «1 

950  END 

960  • 

970  • 

2000  '  LANGMUIR  AST  ALGORITHM 
2100  • 

2110  • 

2115  IF  (PTOT  «  0)  THEN  N1  ■  0:  RETURN 
2120  K1  =  B1  *  EXP (El  /  RT) 

2130  K2  =  B2  *  EXP(E2  /  RT) 

2200  IC  =  0 

2210  PSI  »  .0000001 

2220  FOR  IJ  =  0  TO  1  STEP  1 

2230  PSIG  =  PSI  *  (1  +  DX  *  IJ) 

2240  PIO  *  (EXP(PSIG  /  Ql)  -  1)  /  xi 

2250  P20  =  (EXP (PSIG  /  Q2)  -  1)  /  K2 

2260  XI  *  PI  /  PIO 

2270  X2  =  P2  /  P20 

2280  IF  (IJ  =  0)  THEN  ERl  =  XI  +  X2  -  1  ELSE  DERI  =  XI  +  X2  -  1 
2300  NEXT  IJ 

2310  IF  (ABS(ERl)  <  DX  *  1000)  THEN  GOTO  2500 
2320  IC  »  1C  +  1 

2330  IF  (IC  >  50)  THEN  PRINT  "NON  CONVERGE":  END 
2340  DDERl  =  (DERI  -  ERl)  /  (PSI  *  DX) 

2350  PSI  =  PSI  -  ERl  /  DDERl 
2400  GOTO  2220 
2420  • 

2450  • 

2470  • 

2500  NIO  =  Ql  *  K1  *  PIO  /  (1  +  K1  *  PIO) 

2510  N20  =  Q2  *  K2  *  P20  /  (1  +  K2  *  P20) 

2520  NTOT  =  1  /  (XI  /  NIO  +  X2  /  N20) 

2530  N1  =  NTOT  *  XI 
2550  • 

2600  RETURN 


